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Abstract 

In  this  report  we  describe  two  numerical  techniques  for  the 
computation  of  the  boundary  value  of  time-dependent  potential  flows 
in  a  bounded  domain  where  part  of  the  boundary  is  a  free  surface. 
The  linearized  free  surface  condition  relates  the  normal  derivative 
of  the  potential  to  time  derivatives  of  the  potential  on  the 
undisturbed  free  surface.   It  is  assumed  that  on  the  fixed  part 
of  the  boundary  the  normal  derivative  of  the  potential  is  a  given 
function  of  time. 

The  problem  is  formulated  via  Green's  identity.   The  first 
technique,  QUAKE,  uses  the  classical  time-dependent  Green's 
function  which  satisfies  the  linearized  free  surface  condition, 
and  the  boundary  value  of  the  potential  is  obtained  as  a  solution 
of  an  integral  equation.   The  second  technique,  TDIET,  uses  the 
time-independent  fundamental  solution  of  the  Laplace  equation  and 
the  boimdary  value  of  the  potential  is  obtained  as  a  solution  to 
a  differentio- integral  equation. 
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1.  Introduction 

In  this  report  we  consider  the  following  two-dimensional 
potential  flow  problem:   One  or  several  platforms  are  floating  in 
a  channel;  the  platforms  and  the  channel  are  assumed  to  be 
infinitely  long  in  the  y-direction  and  to  have  a  uniform  and 
bounded  x-z  cross-section  :0  ,   The  fluid  is  initially  at  rest  and 
its  motion  is  generated  by  the  excitation  of  the  confining  boundary 
or  the  platforms.   We  assume  both  the  motions  of  the  fluid  and  the 
platforms  to  be  small  and  we  approximate  the  problem  by  using  the 
two-dimensional  transient  linear  wave  theory. 

Without  loss  of  generality  we  assume  the  boundary  r  of  i^  to 
be  composed  of  coordinate  lines  with  the  undisturbed  water  level 
at  z  =  0;  the  flat  bottom  is  located  at  z  =  -H  (see  Fig.  1). 

The  motion  of  the  fluid,  which  is  assumed  irrotational,  may 
be  described  by  means  of  a  velocity  potential  ^{x,t),    x  =  (x, z), 
satisfying  Laplace's  equation  (see  [4],  [5]) 

(1.1)  ^'^  =   Kx'^Kz   ^  ^   '  ^^  O  >         t  >  0 

and  the  linearized  boundary  condition  on  the  undisturbed  water 
surface  rrn  (z  =  O) 

r 

(1.2)  (l)^^(x,t)  +g(l)^  =  0,  L^   ^Y   '  t>0. 

If  the  free  surface  is  described  by  z  =  Z(x, t),  then 

(1.3)  Z(x,t)  =  -  ^  ^^{x,Q,t)    . 


On  Ta  =  r  -Fp  we  require  the  normal  velocity  -^  (f)(x,  t)  of 
the  fluid  to  coincide  with  a  given  normal  velocity  of  the  rigid 
boundary,  i.e. 

(1.4)  -^  (l)(x,t)  =  a^(x)T(t)  ,  ^  ^  ^k  '         "t  >  0  . 

Here  n  =  (n  ,n  )  is  the  unit  normal  pointing  into  the  domain.   The 
fluid  is  assumed  to  start  from  rest 

(1.5)  Z(x,0)  =  Z^(x,0)  =  0  ,  ■  X  ^  Tp  . 

The  pressure  P(x, z,  t)  is  determined  by  a  modified  form  of 
Euler's  integral 

(1.6)  P(x,t)  +p(f)^  +pgz  =  0  . 

The  term  -^   p|v<{)|  ,  being  of  higher  order,  has  been  sup- 
pressed; the  term  p(x,  t)  =  -pk-\-   is  the  hydrodynamic  pressure. 

Our  problem,  defined  by  (l.l)  -(1.6)  can  be  solved  directly 
by  using  a  finite-difference  approximation  to  the  Laplacian 
operator,  (l.l)  and  the  boundary  condition,  (1.3)  (see  [1]).   How- 
ever in  many  engineering  problems  only  the  solution  on  the 
boundary,  namely  the  wave-height  (1.3)  a-s  well  as  the  pressure 
exerted  on  the  floating  platforms  and  the  confining  walls,  is  of 
interest.   Thus  only  the  solution  of  the  boundary  potential  need 
be  obtained;  a  method  based  on  Green' s  identity  is  tailored  for 
this  situation.   Sometimes  it  is  convenient  to  solve  directly  for 
the  acceleration  potential  ^  =  (j),  .   In  terms  of  the  acceleration 
potential,  the  free  surface  elevation  (1.3)  is  given  by 


(1.7)  Z(x,t)  =  -  I  ^(x,0,t) 
and  the  hydro dynamic  pressure  is  given  by 

(1.8)  P(x,t)  =  -pV/(x,t)  . 

It  is  easy  to  verify  that  ij/  satisfies  Laplace's  equation  and 
the  following  boundary  conditions  which  are  analogous  to  (1.2)  and 
(1.^): 

(1.9)  ^^^(x,t)  +gilJ^{x,t)    =   0   ,  L  ^   ^F   '         t>0, 

(1.10)  -^  ^(x,t)  =  a^(x)T(t)   ,  ^^^A'         ^1^ 
where 

T(t)  =^  T(t)  . 

We  shall  now  describe  two  numerical  techniques  for  the 
computation  of  the  boundary  potential.   In  Section  2  we  shall 
describe  QUAKE,  which  is  a  technique  based  on  the  particular  choice 
of  the  classical  time-dependent  Green's  function;  this  choice 
results  in  an  integral  equation  for  the  boundary  potential,   QUAKE 
is  applicable  to  open  sea  problems  (unbounded  domain)  as  well.   In 
Section  3  we  shall  describe  TDIET  which  is  based  on  a  particular 
choice  of  a  Green's  function  which  is  independent  of  time;  this 
choice  results  in  an  integro-differential  equation  for  the  boundary 
potential.   In  its  present  form  TDIET  is  applicable  only  to  problems 
in  a  bounded  domain.   In  this  case  it  is  much  more  efficient  than 
QUAKE  and  the  finite-difference  methods  in  [1]. 


Both  QUAKE  and  TDIET  are  easily  extendible  to  domains  with 
curved  boundaries  and  to  problems  with  a  three-dimensional  geometry. 
In  Section  4  we  present  some  numerical  results. 


(1.7)  Z(x,t)  =  -  ^  ^(x,0,t) 

o 

and  the  hydrodynajnic  pressure  is  given  by 

(1.8)  p(x.-t)  =  -p^(x,t)  . 

It  is  easy  to  verify  that  ip  satisfies  Laplace's  equation  and 
the  following  boundary  conditions  which  are  analogous  to  (1.2)  and 
(1.4): 

(1.9)  ^^^(x,t)  +giJ^{x,t)    =   0   ,  x^Tp,    t>0, 

(1.10)  -^  ^(x,t)  =  a^(x)T(t)   ,  ^^^A'         "^  -  ° 
where 

T(t)  =^  T(t)  . 

We  shall  now  describe  two  numerical  techniques  for  the 
computation  of  the  boundary  potential.   In  Section  2  we  shall 
describe  QUAKE,  which  is  a  technique  based  on  the  particular  choice 
of  the  classical  time-dependent  Green's  function;  this  choice 
results  in  an  integral  equation  for  the  boundary  potential.   QUAKE 
is  applicable  to  open  sea  problems  (unbounded  domain)  as  well.   In 
Section  3  we  shall  describe  TDIET  which  is  based  on  a  particular 
choice  of  a  Green's  function  which  is  independent  of  time;  this 
choice  results  in  an  integro-dif ferential  equation  for  the  boundary 
potential.   In  its  present  form  TDIET  is  applicable  only  to  problems 
in  a  bounded  domain.   In  this  case  it  is  much  more  efficient  than 
QUAKE  and  the  finite-difference  methods  in  [1]. 
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Both  QUAKE  and  TDIET  are  easily  extendible  to  domains  with 
curved  boundaries  and  to  problems  with  a  three-dimensional  geometry, 
In  Section  4  we  present  some  numerical  results. 


2.    Time-dependent  Green's   F\mction  Method    (QUAKE) 

Let  G(x,_|_,  t),  X  =  (x,  z),  I  =  ii,C.),    be  a  time-dependent 
Green's  function  which  satisfies  the  following: 


(2.1a) 


G{x,i,t)   =  log  R  +W(x,|,t) 


where 


r2  =  (x-e)^+(z-C)^  , 


(2.1b)    A^W  =  0  , 


for  X  ^  <C'  ^^^   t  >  0 


(2.2) 


^tt  ■^g'^z  =  °  ' 


for  ^  e  Tth  and  t  >  0 
—     r 


(2.3) 


G^  =  0  , 


at  ;;  =  -H  and  for  t  >  0 


(2.4a)    G(x,^,t)  =  G(x,e,-t)  , 


for  all  t 


(2.4b)    G(x,^,0)  =  0  , 


for  ^  e  r. 


Such  G(x,_|_,  t)  exists  and  is  given  explicitly  by  the  following 
expression: 


oo 


fP   p:a^      r   -   inp    fR/Rn      P    /      r^"   sinh  (Kz)sinh  (kc)cos   (Kr)  dK 
(2.5a)      G  -    log    (R/R    j  -2  J      K  cosh(KH)  

0 

oo 
.    2    f       cosh  K(z+H)cosh  K(C+H)  [1  -cos    (u3t)]     ^^^   (Kr)dK 
^  K   cosh''^(KH)   tanh(KH) 

where 


2.2 


(2.5b)      cn"^    =   gK  tanh(KH)    ,       R'    =    [  (x-^""  +(z+C)    1       .       ^   =    |x-||     . 


(See  [2],    [4]  and  [5].)   Observe  that,  as  one  expects,  G(x,^, t)  is 
symmetric  in  x  and  |. 


In  order  to  derive  an  integral  equation,  we  apply  Green's 
theorem  to  the  functions  ^(x,t)  and  G(x,^,t-T)  (see  [5]) 

(2.6)  27r^(x,t)  =^  [G(x,i,t-T)^^(i,t)  -G^(x,e,t-T)^(4,t)]ds^ 

r 

for  all  X  e  jj 
Integrating  (2.6)  with  respect  to  t  between  t  =  0  and  t  =  t  we  find 

T 

(2.7)  27r[<t)(x,T)  -<|)(x,0)]  =J   j)   [G(x,i,t-T)^^(i,t) 


0  r 

-  G^(x,e,t-T)^(e,t)]ds^dt  . 

We  observe  that  because  of  the  boundary  conditions  (1.9)  and 

(2.2),  the  integrand  on  r-o  in  (2.7)  can  be  written  as  the  following 

r 

time  derivative 

Using  this  relation  together  with  (1.5 )  and  (2.4),  we  get 

T 

(2.8)  J   j^  [G(x,i,t-T)^^(i,  t)  -G^(x,e,t-T)^(e,t)]ds^ 
1    ^ 


Ij  [[G(x,i,0)^^(e,T)  -G^(x,|_,0)^(i,T)] 

[G(x,|_,T)^^(e,0)  -G^(x,^,T)^(e,0)]]ds^  =  0 


^F 


Thus 


the  integral  in  (2.7)  can  be  only  taken  along  Ta  =  r  -r-p» 


Substituting  for  ijj      on  Ta  the  values  (l.lO),  we  obtain 
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(2.9)  27r[Mx,T)  -Mx,0)]  =J   j^  [G(x,i,t-T)a^(|_)T(t) 

Now  we  differentiate  equation  (2.9)  with  respect  to  x  and 

use  the  relations  G^   =  -G^,  G^^.  =  -Gj^^  and  (2.4) 

(2.10)  2TrV'(x,T)  =y^  [G(x,e,0)a^(_|_)T(T)  -G^(x,i,0)^(|_,T)]ds^ 

+  /  /  [G^^(x,e,t-T)^(l_,t) 

-   G^(x,i,  t-T)a^(_|)T(t)]ds^dt  . 

Next  we  let  x  e  JJ      in  (2.10)  approach  a  point  on  the  boundary  Ta. 
This  results  in  the  following  integral  equation  for  the  boundary- 
potential: 


(2,11a)   a(x)^(x,t)  +  P.V.  /  G^(x,_|_,  0)^(_|_,  T)ds 


I 


Ta 


=  J     J   Gj^^(x,|_,  t-T)^(|_,  t)ds^dt  +  c(x,t)  , 


0  r^ 


(2.11b)   c(x,t)  =  -J    G(x,e,0)a^(^)T(T)ds^ 


Ta 


T 

+  j^  [T(t)  J  G^(x,|_,t-T)a^(^)ds^]dt 


°       ^A 

Here  P.V.  indicates  the  Cauchy  principal  value  integral  and  A(x)/27r 
is  the  part  of  an  infinitesimal  circle  centered  at  the  boundary 
point  x,  which  is  contained  in  .^)  . 


The  convolution  with  respect  to  time  in  the  right-hand  side 
of  (2.11)  is  approximated  by  some  quadrature  formula,  e.g.  the 
trapizoidal  rule: 

T 


0  r 


A 


r 


A 


L-l 


+  At  ^ 


&=1 


^nt^ii'i'-^^^^^^i'^^-^)^*)^^! 


where  L  =  x/At. 

Next  we  discretize  (2.11)  with  respect  to  space  by  dividing 
the  boundary  Ta  into  intervals  A.,  1  Jf  i  J^  N,.   QUAKE  uses  a 
piecewise  constant  approximation  to  ■>{/;    the  constant  values  are 
■>(/.    =  Tp{x. ,  £At)   where  x.  is  the  center  of  A.. 

With  this  discretization  (2.11)  becomes  a  system  of  N  linear 
equations 


(2.12) 


where 


A^^  =  At  XZ  bV^'^  +7^ 
to   ~ 


\^i 


,^N 


an.d  A  is  an  N„  X  N.  constant  matrix  whose  entries  are 
(2.13a)      A.  .  =  5.  .a(x.  )  +  P.V.  /  G  (x.,^,0)ds. 


A, 


At 
2 


B   is  an  N»  xN.  matrix  whose  entries  are 


^  is  an  N,  vector  whose  components  are 

(2.l4a)  -y^  =   c(x^,lAt)    =   -  a^T(rAt)  -^[T(0)p[' +T(lAt(p°] 

L-1 

-   At  Y~  T((L-£)At)pf 

i^  ^ 

where  we   have   introduced 

(2ol4b)  a^   =jrG(x^,e,0)a^(e)ds^    , 

(2.14c)  P^   =  J"G^(x^,e,  Mt)a^(Ods^    . 

^A 

We  observe  that  the  matrix  A  and  the  vector  a  are  independent 
of  time  and  therefore  have  to  be  computed  only  once.   The  number  of 

operations  needed  to  compute  L  time-steps  of  the  potential  on  Ta  is 

2   3 
proportional  to  L  N. .   The  quadratic  dependence  on  the  number  of 

time  steps  is  due  to  the  need  to  evaluate  the  convolution  integral 

in  (2.11).   Obviously  QUAKE  requires  a  large  memory  allocation 

because  of  this  convolution.   If  the  free  surface  potential  is 

needed,  then  an  additional  computing  time  has  to  be  expected. 

QUAKE  can  handle  unbounded  domains  as  well  as  bounded 

domains.   In  fact,  it  is  less  costly  for  unbounded  domains  than 

for  a  corresponding  problem  in  a  bounded  domain. 


3.  TDIET  (Time-dependent  Dif ferentlo-Integral  Equation  Technique) 

Let  G(x,|_)  be  the  following  time-independent  Green's  function 

(3.1)  G(x,i)  =  log  R  +W(x,i) 

where  R  is  defined  in  (2.5b)  and  W(x,|_)  is  any  function  which  is 
harmonic  for  all  x  e  J^ .      ^(x, t)  satisfies  Laplace's  equation  in 
L    for  all  t  >  0.   Therefore,  the  following  Green's  identity  is 
satisfied  for  all  t  >  0 

(3.2)  A(x)  ^(x,t)  =(j>   [G(x,i)  -^  '//(i,t)  -^(i,t)  ^  G(x,_|)]ds^ 

r 

where  A(x)/2Tr  is  the  part  of  an  infinitesimal  circle  centered  at  x 

which  is  contained  in  -O  ,      Substituting  the  boundary  conditions 

(1.9)  and  (I.IO)  for  -2^  in  (3o2)  results  in  the  following  integro- 

dn 

differential  equation 

(3.3)  A(x)^(x,t)  +(^^(_|,t)  -^  G(x,e)cis^ 


+  "l/  G(x,i)^^^(i,t)ds^  =  T(t)  J  G(x,e)a^(i)ds^ 


A  convenient  choice  of  a  Green's  function  is  one  for  which 


■5H 


G(x,_|_)  =  0  for  _!_  e  Fa'   In  this  case,  the  second  integral  on  the 


left-hand  side  of  (3.3)  is  taken  only  along  Yt^   and  the  resulting 

discretization  leads  to  a  linear  system  of  ordinary  differential 

equations 

,2 
(3.4)  A  ^  ^p(t)  +B  V^(t)  =  a(t) 

dt 
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for  the  Np  discrete  values  of  the  potential  ^  along  the  free 
surface  r^;  A  and  B  are  N„  x N„  matrices.   After  obtaining  an 
approximatic:     ^,  t)  for  the  potential  on  r^  hy  solving  (3.4)  as 
an  initial  value  problem,  the  approximate  solution  ^. (x, t)  for  the 
potential  on  Ta  =  r  -Tp  is  obtained  from  solving  (3.3)  with 
G{x,i)   =  log  R, 

(3.5)  A(x)V^^(x,t)  +  j  ?^(i,t)  ^  log  R  ds^ 

=  T(t)  Jlog  R.a^(e)ds^  -/[^F^i'"^^  Si  ^°^  ^ 

1         ^2 
+  I  log  R  ._2_  ^(|,t)]ds   . 
S        at^  ^  -      ^ 

The  right-hand  side  of  (3.5)  is  known  and  therefore  the  discretiza- 
tion of  (3.5)  leads  to  the  following  system  of  linear  equations 

(3.6)  C  '±^{t)  =   b(t) 

for  the  N.  discrete  values  of  the  potential  along  r ^i    the  N ,  x N . 
matrix,  C,  is  independent  of  time. 

This  approach  cannot  be  implemented  directly  in  a  numerical 
code  since  an  explicit  expression  for  G(x,_4)  such  that 
-jS.  g(x,  ^)  =  0  for  ^  e  Ta^  is  not  always  available.   In  the 
following  we  shall  present  a  technique  to  arrive  at  (3. 4)  and  (3*6) 
without  requiring  an  explicit  form  of  the  Green's  function.   The 
main  idea  in  TDIET  is  to  first  discretize  (3*3)  with  G(x,J_)  =log  R 
and  then  to  obtain  the  desired  decomposition  by  an  ordered 
elimination  procedure. 
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(i)  Discretization:   Let  x.,  1  ^  J  Jf  N,  be  a  partition  of  the 

boundary  and  let  ^  .  (t )  =  ^(x  ., t ) .   It  is  convenient  to  order  the 

J         J 

partition  points  so  that  x.,  1  <  j  <  N»,  are  the  points  along  Tn 

and  X.,  N.+l  jf  J  J^  N,  are  points  along  the  free  surface  T-pi 

N„  =  N  -N..   We  denote  by  S.(x)  the  unit  cubic  spline  such  that 
■f        A  J  — 

(3.7)    Sj(x^)  =  5.  J     (5^j  =  0  for  i  ?^  j,  b^.   =  1  for  i  =  j)  . 

In  the  particular  geometry  of  Figure  1,  the  boundary  has  corners. 
We  construct  S.(x)  so  that  it  will  be  continuous  up  to  the  second 
derivative  in  a  linear  segment  which  includes  x.  and  to  be 
identically  zero  in  all  other  segments.   This  way,  the  approxima- 
tions 

N 
(3.8a)  f(x,t)  ^^S  (x)^.(t)  , 

N  ,2 

(3.8b)  ^   (x,t)  ^  >—  S  (x)  -^  ^  (t) 

^^  ~      j^   J  ~   dt^   J 

are  continuous  up  to  the  second  derivative,  except  at  the  corner 
points  where  only  continuity  of  the  function  is  assumed.   Because 
of  the  singular  nature  of  the  corners,  the  partition  points  are  so 
distributed  that  Iax.I  =  |x.,,  -x.l  is  smaller  in  the  neighborhood 
of  the  corners. 

Next,  equation  (3.3)  is  discretized  by  using  the  cubic  spline 
approximation  (3.8) 


12 


N      r     ^ 

(3.9)    A(x^)^^(t)  +  ^  i^.it)  J   S.ii)   ^  G(x^,e)ci£ 


r 

N     .2 


i  HZ  ^^.(t)  /  S.(OG(x.,e)d£ 


=  T{t)  J  G{x^,i)a^{i)ds^    , 


1  <  i  <  N  . 


The  N  equations  (3.9)  can  be  written  in  the  following  matrix  form 


(3.10) 


E.M.(t)  =  T(t)a 


where  p.  is  an  N  +N„  vector  whose  components  are 


(3.11) 


i^i(t)  = 


f  ■^A^-i'"'^^     ' 

2 

^^(x^,t)     , 


ht' 


^F^^i-Np'^)  ' 


1  <  i  <  N^ 


N^+1  <  i  <  N 


N+1  <  i  <  N  +  N, 


a(t)  is  an  N  vector  the  components  of  which  are 


(3.12)     a^  =  jG(x.,|)a^(i)ds 


i 


1  <  i  <  N 


E  is  an  Nx(N+N-p)  matrix  whose  entries  are 

(3.13a)    E^^.  =  A(x^)5^j  +/Sj(i)  -^  G(x^,e)ds^  ,    1  ^  1  ^  N  , 

r 

(3.13b)    E^^.  =1/  Sj(|_)G(x^,|)ds^  ,    1  ^  i  ^  N  ,   N^+1  <  j  <  N  , 
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p  • 


(J.Uc)   E.J  =  ^(il)6l^J.„  +/Sj_„  (1)  ^  G(x.,l)as^  , 

r 

1  <  i  <  N  ,   W+1  <  ,]•  <  N  +N. 

We  ol:       that  E.  .,  1  <  i  <  N,  1  <  j  <  N  +N„,  as  well  as  a., 
1  <  i  <   .  .  pend  only  on  the  geometry  of  the  boundary  and  are 
indep^,.^.,..t  of  time,  provided  that  the  Green's  function  is  inde- 

:'  timeo      ,  ,   •  .     ■  -   '■    earlier,  we  choose 
G(x,  O  =  1-og  Rj  R^  =  (x-^)"^  +(z-^)'^,  for  an  arbitrarily  shaped 
domain.   In  case  the  bottom  is  flat  as  in  Figure  1,  we  take 
G(x,|,)  =  log  RR,  R^  =  (x-e)^  +(z+C  +2H)^.   This  way  -^  G(x,|)  =  0 


on  the  flat  bottom  C,   =  -H.   Therefore  the  integrand  in  the  integrals 
in  (3.13)  is  identically  zero  on  the  bottom,  and  the  potential 
values  on  the  bottom  are  eliminated  from  the  equations  (3.10)- 

(3  .  13  ) . 

We  remark  that  in  the  particular  geometry  of  :  -  .  :  1,  where 
the  boundary  is  composed  of  coordinate  lines,  the  coefficients 
(3.12)- (3.13 )  may  be  computed  analytically.   This  is  needed,  at 
least  locally,  because  of  the  singularity  of  the  integral 
x^  =  !_  (see  [3]  )o 

(ii)  Solution  of  the  discrete  problem:   In  order  to  obtain  an 

an-:-  ,^  ;_   ::e  equations  (3.4)  and  (3.6),  a  Gauss  elimination  (with 

row  T  ■   '  '  ■  :)  is  ■  :"  ■ -^d  to  the  matrix  equation  (3. 10):   The 

w  ;,...,u.  -I  are  e.^  '       1  from  the  j-th  equation, 

2  <  j  <  N.   X    -      ■  _  '  n  is  of  the  following  form 
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<-  N^  ^   <-  Np  -   <-  Np  ^ 


(3.14) 


t 
N, 


N. 


"a  i 

1 
■ 
I 
• 

Q 

1 
1 

R 

0 
+  _ 

""u; 

_4_. 

1 

0      : 

1 
1 

S 

1 

0 

n' 

l^it) 


ip(t) 


l^it) 


-A 


^F 


T(t) 


The  last  N„  equations 


(3.15) 


"f  "^  S^^)  +SV^(t)  =  c^  hT  T(t) 


-F  dl 


comprise  a  sy 


f  linear  ordinary  differential  equations  similar 


to  (3.4);  U„  is  an  upper  diagonal  matrix.   This  subproblem  for  the 
r 

discrete  values  of  the  potential  on  the  free  surface  is  now  solved 
independently  as  an  initial  value  problem  in  the  time  interval  of 
interest.   The  initial  value  problems  (3.15)  and  (1.5)  can  be 
solved  numerically  by  any  standard  method;  our  computer  code  uses 
a  fourth  order  accurate  method  of  Runge-Kutta. 

Once  the  free  surface  potential  *-p(t)  has  been  computed  in 
the  time  interval  of  interest,  '^A('t)  can  be  found  at  the  desired 
time-levels  by  solving  the  system  composed  of  the  first  N.  linear 
equations  in  (3.l4) 


(3.16) 


U^3^^(t)  =  ^T(t)c^  -Q  ^  ^p(t)  -R^p(t) 

LI  U 


Observe  that  U.  is  an  upper  diagonal  matrix.   Hence,  once  the 
right-hand  side  of  (3.l6)  is  computed,  the  solution  i//„(t)  is 
obtained  by  a  simple  back  substitution. 
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In  its  present  form,  TDIET  is  applicable  for  problems  in  a 
bounded  domain  where  only  the  boundary  potential  is  needed.   In 
case  the  potential  is  needed  for  a  few  time  levels  in  the  domain, 
it  still  is  efficient  to  use  TDIET  to  compute  the  boundary 
potential  at  the  desired  time  levels  and  then  to  solve  a  Dirichlet 
ptoblem. 

We  have  described  TDIET  for  two  dimensional  computations. 
For  the  formulation  in  three  dimensions  one  has  only  to  change  the 
Green's  function  and  the  notation. 
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4.  Numerical  results 

In  this  section  we  shall  present  several  numerical  results 
which  will  demonstrate  the  efficiency  and  the  accuracy  of  TDIET. 

(i )  Free  surface  disturbances 

We  consider  an  empty  rectangular  basin  0  <  x  <  2L, 

-H  _<  z  _<  0,  in  which  the  confining  walls  and  the  bottom  are  at 
rest. 

(4.1)  (t)^(0,z,t)  =  i^(2L,z,t)  =(t)^(x,-H,t)  =  0  ,    t  >  0  . 

The  motion  of  the  fluid  is  generated  by  initially  disturbing  the 
free  surface 

00 

(4.2a)  <t>(x,0,0)  =  I(x)  =}         I^  cos  iJg  X  , 

n=0 

OP 

(4.2b)  (J).  (x,0,0)  =  J(x)  =  y^   J  cos  ^  X 

or  by  applying  a  periodic  pressure  impulse  to  the  free  surface 

°°  n-rr 

(4.2c)  P(x,  t)   =  p(x)    sin  Ot  =   sin  Ot  ^         P^   cos  -IL  x    . 

The  free  surface  condition  (1.6)  becomes 

(4.3)  i^^  +g<^^   "  "  ?  ^t^^'^^  "  "  P  ^^^''  ^°^  "^ 

and  a  corresponding  inhomogeneous  term  is  added  to  the  system  of 
the  ordinary  differential  equations  (3. 10). 

It  is  easily  verified  that  the  solution  to  the  problem  (4.1)- 

(4.2)  is 
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00     cosh  ^  (z+H) 

(4.4a)   ())(x,z,t)  =  }        I cos  ^  x  cos  o)  t 

T^     1^   cosh  ^  H       ^ 


n 


+  J  t  +  5    J 


00     cosh  ^    (z+H) 


n=l    o)  cosh  -r-  li 
n       L 


^   00     cosot-  coso)  t  cosh-r=— (z+H) 

+  Qy—  P  , — __i} LJ Icosl^x, 


n 


cosh  ^  H 


(4.4b)   Z(x,t) 


J     -,   oo 

—  +  —  J         cos  -^  X  (l  (D   sin  CO  t  -  J   cos  co  t 
g    g  ^-^y      L    ^  n  n      n    n      n 


where 


^    O  sin  Ot  -GL!   sin  (x;  t 
,  O  p  n r^N 

p   n        ~2    2       ^ 


(4.4c) 


a   =  -T-  ,    (X:   =  ga   tanh  a^H 
n    L  '     n   '^  n       n 


Example  1:   Initial  evaluation 


(4.5a) 
(4.5b) 


(t)(x,0,0)  =  0  , 


49 


^^(x,0,0)  =    lOYZ   ^  cos  [(2n+l)7r  ^^^] 


n=0  (2n+l)' 
'   5u  ^^+50)  , 


J^(50-x)  , 


■  100  <■  X  <  0 


0  <  X  <  100   , 


(4.5c)     P(x,t)  =    0    . 

The  initial  free  surface  elevation  is  shovm  in  Fig.  2a;  25 
points  T:        :  ly  distributed  along  the  undisturbed  free  surface 
z  =  0,  -100  <  X  <  100.   The  depth  is  H  =  10.   This  problem  is 
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<   t    <   20,  :ne    steps.      Figures  2a,    2t),    2c, 

2_    _____    ^                  _:_-:    calcula_    -  ...                         n  at  t  =  0,    8, 

12,    l6,    2r~     -           '  ■     -"     „      Tatle  1    r'     v:^    a  comparison  b'-^-..-' -n   the 

com].  1    (4.4).  -elative 


L^-  t   =    c 

c 


25  '     25 


25         .„^    "     .._^    o/     25  T-^ 


(4.6a)      e(^)   =    j  J^   [^™-i^"-       .   ^'   J^   [^i'         .   j 

(4.6b)      e(Z)    =    [;ri    [Z^^-Z^^^^h^/^    ^2f^ct^2l2  ^   ^^^^ 


The  CP  time  on  the  CDC  6600  for  this  computation  is  5*85  seconds: 
2.32  seconds  to  arrive  at  the  discrete  formulation  and  1.53  seconds 
to  advance  the  solution  in  50  time  steps. 

Example  2:  -ure  ir  • 

(4.7a)   (j)(x,0,0)  =  (|)^(x,0,0)  =  0  ,  -100  <  x  <  100 

(4.7b)   P(x,  t)  =  sin  0.47rt  .p(x)  ,  -100  <  x  <  100 

(4.7c)   p(x)  =  g  +  4\ZI  [(4-"^^^)  sin(naTr) 

-^   T  a  n=l   n 

_  2HL  cos  (naTr)] 
n      \    ^  J 


^^lOOa^    -^J   ' 


0 


(-1)'^ 

cos     (nTT 

X  +100 

n5 

100 

X 

<   100a 

X 

>   100a 

p(x)  given  by  (4o7c)  with  a  =  0.25  is  shown  in  Figure  Ja.  The 
geometry  in  this  example  is  identical  to  the  one  in  Example  1, 
This  problem  is  solved  for  0  <  t  <  10  (2  periods),  using  80  time 
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steps.   Figures  2t),  2c,  2d  and  2e  show  the  calculated  free  surface 
elevation  (the.  continuous  line)  and  the  exact  free  surface  eleva- 
tion (4.4b)  (asterisks)  at  t  =  2.5,  5.,  7.5^  10.,  respectively. 
The  relative  Lp-errors  (4,6)  at  t  =  10  are  e((j))  =  0.425^  and 
e(z)  =  0.6145^.   The  CP  time  on  the  CDC  660O  for  this  calculation 
is  4.8  seconds:   2.4  seconds  to  arrive  at  the  discrete  formulation 
and  2.4  seconds  to  advance  the  solution  in  80  time  steps. 

(ii)  Harmonic  excitation  of  walls 

We  consider  an  empty  rectangular  basin  0  <  x  _<  2L, 
-H  <  z  <  0,  in  which  the  bottom  is  at  rest  and  the  free  surface  is 
undisturbed. 

(4.6a)     (|)(x,0,0)  =  (f),  (x,0,0)  =  0  ,    0  <  X  <  2L  , 

(4.6b)     P(x,t)  =  (|)^(x,-H,t)  =0   ,    0<x<2L,    t>0. 

The  motion  of  the  fluid  is  generated  by  applying  harmonic 
excitation  to  the  walls 

(4.7)   (f)  (0,  z,t)  =  (j)  (2L,z,t)  =  -A  cosOt  ,   -H<z<0,   t>0. 

A.  tC  —     —  — 

The  solution  to  this  problem  can  be  approximated  by  the 
following  expression 


4a         00    P 
(4.8a)   tt>(x,z,t)  =  -^  cos  Ot  YH  ^Iti+l 


1  - 


^2    cosh  ct^]^_^_^{z+E)' 

TJ^        2      cosh  a„,  ,-,H 
O  -0)2^,+  !         2k+l   . 


cos  (cx2j^_^^x) 


p 
4A  ^2-   _2      ^Vk+1    ^°^^   °'2k+l^^+^) 
+  ^  IZ  -2^+1  -^^ cosh  SI..  H    ^°^  «2k+l-  ^°^  -n^ 


k=0  "^^  -^  O  -0)21^+1        ^2k+r 
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where  ou     and  ax     are  given  by  (4.4b).   The  free  surface  elevation 
is  given  by  the  following  expression 


_2 
1   - 


7^ — 2~ 

O  -">2k+l-J 


cos  ^2]!i+'!L^ 


4a  CO   _p 

(4.8b)   Z(x,t)  =  g  n  sin  Ot  y^  g^g^^ 

4a   '-^   -2     "^k+1 

+  gl  ^  °'2k+l  72 5 ^°"  °'2k+l^  ^i^  ^2k+l^   . 

k=0       O  -CD2i,+i 

Figure  4  shows  a  comparison  between  (4.8b)  and  the  solution 
of  TDIET  with  o  =  2Tr  at  t  =  0.75.   The  geometry  of  the  basin  is 
given  in  the  figure. 

(iii)  Heaving  of  a  rectangular  cylinder 

In  this  example  we  compute  the  forced  heaving  motion  of  a 
semi-submerged  rectangular  cylinder  in  the  free  surface.   As  in 
[1]  we  choose  a  rectangular  cylinder  with  a  draft  D  =  1  and  a  beam 
B  =  2;  the  water  depth  is  H  =  7.15  (see  Fig.  l).   The  vertical 
velocity  of  the  cylinder  is 

(4,9)        i  (x, -D,  t)  =  O.IO  cos  Ot  ,    O  =  1.253  /g  . 

In  order  to  simulate  open  sea  conditions,  the  walls  are  placed  at 
X  =  ± l6'  Although  the  solution  is  symmetric  in  x  we  solve  the 
complete  problem.   Each  of  the  free  surface  segments  1  _<  |x|  _^  l6 
is  divided  into  24  unequal  intervals,  the  size  of  which  increases 
quadratically  with  |x|.   The  bottom  of  the  cylinder  is  divided 
into  12  intervals.   The  sides  of  the  cylinder  and  the  walls  are 
each  divided  into  8  intervals. 

Figures  5a,  5b  and  5c  show  the  wave  profile  at  the  end  of 
1,  2,  and  three  oscillations  of  the  cylinder,  respectively. 
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Observe  the  establishment  of  a  standing  wave  pattern  near  the  body 
and  that  the  asymptotic  wave  length  A   =  27rg/o  =  4  is  almost 
reached.   These  wave  profiles,  although  qualitatively  similar,  are 
not  identical  to  those  shown  in  [1];  this  is  probably  due  to 
different  initialization  of  the  problem.   However,  the  dynamic 
pressure  force  agrees  quite  well  with  the  results  reported  in  [1]. 

The  CP  time  on  the  CDC  6600  for  this  calculation  is  23 
seconds:   12  seconds  to  set  up  the  discrete  formulation  and  11 
seconds  to  advance  the  solution  in  90  time  steps  (30  time  steps  per 
oscillation) . 

(iv )  Seismic  excitation  of  a  rectangular  cylinder 

In  this  last  example  we  compute  the  forces  exerted  on  a  semi- 
submerged  rectangular  cylinder  which  is  held  fixed  in  the  fluid. 
The  horizontal  confining  walls  are  excited  by  an  accelera- 
tion which  resembles  earthquaks  data  (Fig.  6a) «   This  problem  is 
solved  by  both  TDIET  and  QUAKE  and  the  time  history  of  the  horizon- 
tal hydrodynamic  force  on  the  cylinder  is  compared  (Fig.  6b;  the 
geometry  is  shown  in  the  figure).   Both  methods  respond  in  a  very 
similar  way  to  the  seismic  acceleration.   The  crude  approximation 
used  to  derive  QUAKE  is  probably  responsible  for  the  numerical 
discrepancy  between  the  two  results. 

As  was  earlier  remarked,  the  convolution  operator  in  QUAKE 
makes  it  extremely  expensive  for  long  time  calculations,  e.g. 
a  computation  of  hydrodynamic  forces  due  to  an  earthquake  with  a 
duration  of  50  seconds  by  TDIET,  is  about  600  times  faster  than  the 
corresponding  computation  by  QUAKE  (using  the  same  number  of  time 
steps  and  a  similar  spacial  partition). 
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Table  1:   Comparison  between  numerical  solution 
and  exact  solution  in  Exaimple  1, 


I 

'''numerical 

^exact 

numerical 

2 
exact 

1 

-18.951 

-18. 794 

-.0352 

- , 0284 

2 

-17.979 

-17.880 

-.0335 

-,0408 

3 

-15.521 

-15.478 

-  0  0244 

- . 0195 

4 

-12.082 

-11.997 

- . 0126 

-,0116 

5 

-    8.284 

-   8,259 

-.0161 

- . 0194 

6 

-   4.321 

-   4.315 

-.0179 

-.0172 

7 

0.000 

0.000 

.0000 

,0000 

8 

4.321 

4.315 

.0179 

.0172 

9 

8.284 

8.259 

.0161 

.0194 

10 

12.082 

11.997 

.0126 

.0116 

11 

15.521 

15.478 

.0244 

.0195 

12 

17.979 

17.880 

.0335 

.04  08 

13 

18.931 

18. 794 

.0352 

.0284 

14 

17.979 

17. 880 

.0335 

.04  08 

15 

15.521 

15.478 

.0244 

.0195 

16 

12.082 

11.997 

.0126 

,0116 

17 

8.284 

8.259 

.0161 

,0194 

18 

4.321 

4.315 

.0179 

,0172 

19 

0.000 

0.000 

.0000 

,0000 

20 

-   4,321 

-  4,315 

-.0179 

-.0172 

21 

-    8.284 

-  8.259 

-.0161 

- , 0194 

22 

-12.082 

-11.997 

-.0126 

-.0116 

23 

-15.521 

-15.478 

- . 0244 

-.0195 

24 

-17.979 

-17.880 

-.0335 

-,0408 

25 

-18,931 

-18. 794 

-.0352 

-.0284 
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Figure  2a:   Initial  surface  elevation  (Example  l), 
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Figure  2b:   Surface  elevation  at  t  =  8  (Example  l) 
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Figure  2c:   Surface  elevation  at  t  =  12  (Example  l) 
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Figure  2d:   Surface  elevation  at  t  =  l6  (Example  l), 


28 


ir.  ■ 


a 

(Tl- 


^-ptrtnr' 


cc". .. 
>  I 

UJ 


-60-00 


-20-00 

X 


20.00 


60-00 


lt)0.00 


Figure  2e:   Surface  elevation  at  t  =  20  (Example  l) 
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Figure  3a:   Pressure  impulse  (Example  2) 
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Figure  3b: 


Surface    elevation   at   t   =  2.5    (Example    2) 
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Figure  3c 


Surface   elevation  at   t   =  5    (Example   2) 
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Figure  2e:   Surface  elevation  at  t  =  20  (Example  1), 
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Figure  3a:   Pressure  impulse  (Example  2) 
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Figure  3b 


Surface  elevation  at  t  =   2.5  (Example  2) 
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Figure  3c:   Surface  elevation  at  t  =  5  (Example  2). 
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Figure  3d:   Surface  elevation  at  t  =  7,5  (Example  2) 
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Figure  3e:   Surface  elevation  at  t  =  10  (Exajnple  2), 
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FIGURE   4     SURFACE  ELEVATION    AT     t  =  0.75  SEC 
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{^    Figure   5a 


Free   surface   profile   at   the   end   of  one   oscillation 
(Heaving   cylinder). 
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Figure   5t> 


Free  surface  profile  at  the  end  of  two  oscillations 
(Heaving  cylinder). 
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^1   Figure    5c 


Free  surface  profile  at  the  end  of  three  oscillations 
(Heaving  cylinder). 
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FIGURE  6a     HORIZONTAL   SEISMIC  ACCELERATION  AT  WALLS 
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FIGURE  6b     HORIZONTAL  HYDRODYNAMIC  FORCE 

ON  BODY  INDUCED  BY  AN  EARTHQUAKE. 
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This  report  was  prepared  as  an  account  of 
Government  sponsored  work.   Neither  the 
United  States,  nor  the  Administration, 
nor  any  person  acting  on  behalf  of  the 
Administration: 

A.  Makes  any  warranty  or  representation, 
express  or  implied,  with  respect  to  the 
accuracy,  completeness,  or  usefulness  of 
the  information  contained  in  this  report, 
or  that  the  use  of  any  information, 
apparatus,  method,  or  process  disclosed 
In  this  report  may  not  infringe  privately 
owned  rights;  or 

B.  Assumes  any  liabilities  with  respect  to 
the  use  of,  or  for  damages  resulting  from 
the  use  of  any  Information,  apparatus, 
method,  or  process  disclosed  in  this 
report . 

As  used  in  the  above,  "person  acting  on  behalf 
of  the  Administration"  Includes  any  employee 
or  contractor  of  the  Administration,  or 
employee  of  such  contractor,  to  the  extent 
that  such  employee  or  contractor  of  the 
Administration,  or  employee  of  such  contractor 
prepares,  disseminates,  or  provides  access  to, 
any  Information  pursuant  to  his  employment  or 
contract  with  the  Administration,  or  his 
employment  with  such  contractor. 
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a. 


